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SUMMARY 



The use of computers to optimize free parameters of a system has become 
relatively widespread in many areas of engineering. Parameter optimization 
codes have been written for that purpose, and make it possible for a design 
engineer, once he has developed the mathematics of a system, to optimize its 
parameters according to some criteria. But of equal interest to the design 
engineer is the sensitivity of the optimized criteria to departure of the 
parameters away from their optimal value. The purpose of this thesis is to 
show ways in which a parameter optimization code may be augmented to yield 
such sensitivity information. 

A Fletcher and Powell version of Davidon’s variable metric optimization 
search technique was employed to optimize multi-parameter functions. Their 
method is useful in that it computes the inverse Hessian matrix (or matrix of 
second derivatives) , which completely describes the curvature of the function 
at the optimum. Equations were developed so that the sensitivity could be 
expressed in a meaningful output format. This was made possible through the 
use of matrix inversion and eigenvalue analysis subroutines which v;ere ob- 
tained from the scientific subroutine library of the IBM 360 91 and used in 
conjunction with a digital computer code employing the Fletcher and Pov/ell 
technique. Equality constrained optimization problems were also considered 
by employing the penalty factor method proposed by Courant and used by Kelley. 
Equations analogous to the use of Lagrangian Multipliers were used to deter- 
mine the cost of the equality constraint. 

Example problems are offered showing the optimal solutions, sensitivity 
data at the optimum, and interpretation of that data. The well known 
Rosenbrock function was used to exhibit the accuracy of the methods employed. 
A typical engineering problem was solved involving the sensitivity of an 
optimal nuclear rocket engine used to inject a payload onto an interplan- 
etary trajectory. The results indicated that the thermal power of the re- 
actor and the ratio of length to diameter of the core could be varied con- 
siderably from their optimal values with little cost. The power density 
however was relatively fixed for optimal operation. 
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I. INTRODUCTION 

The use of computers to optimize free parameters of a system has be- 
come relatively widespread in many areas of engineering. Parameter opti- 
mization codes have been written for that purpose, and make it possible 
for a design engineer, once he has developed the mathematics of a system, 
to optimize its parameters according to some criteria. But of equal in- 
terest to the design engineer is the sensitivity of the optimized criteria 
to departure of the parameters away from their optimal value. The purpose 
of this thesis is to show ways in which a parameter optimization code may 
be augmented to yield such sensitivity information. The interest in this 
is evident. 

Since the engineer is working with temperatures, pressures, masses, 
etc., the computed optimal solution cannot be implemented exactly. Physical 
parameters are subject to uncontrollable variations and the resulting de- 
parture from the optimum mhy be quite significant. For space flight 
applications, maximum payloads might be of primary interest for one mission 
while a minimum fuel consumption the criterion for the next. Design require- 
ments for a Venus flyby will certainly differ from those of a Mars lander 
or Jupiter probe, and yet many of the systems must be flexible enough to 
be useful for all three. Thus, the optimized solution must also contain 
significant information about departures from the optimum. This thesis 
addresses that problem; i.e. sensitivity analysis at the optimum. 

Parameter optimization algorithms are many in number and varied in 
application. They differ primarily in their rate of convergence and the 
restrictions imposed on the function under consideration. Nearly all 
require a large number of iterations to achieve a given accuracy in locating 



the optimum, and some procedures may not converge from a poor starting point. 

Because, near the optimum, the second order terms in the Taylor series 

expansion dominate, the only method which will converge quickly for a gen- 
eral function are those which will guarantee to find the optimum of a 
general quadratic function speedily. Fletcher and Powell [1]*'^ have pro- 
duced such a method which was based upon a procedure described by Davidon 
[2]. The method is superior to other techniques which possess quadratic 
convergence in that it makes use of information determined by previous 
iterations and also in that each iteration is quick and simple to carry 
out. The primary justification for its use however, lies in the fact that 
the method yields the necessary information to determine the curvature of 
the objective function at the optimum. 

The method of Fletcher and Powell estimates and continually updates 
the inverse of the Hessian matrix, H, (or matrix of second derivatives) 
during the optimization search so that a close approximation to the true 
value at the optimum is reached. Through the use of the eigenvalues and 
eigenvectors of the Hessian matrix, analysis of the characteristics of the 
objective function in the neighborhood of the optimum is possible. 

As previously mentioned, it is of interest to know not only where the 
optimum is located, but by how much each parameter may be changed from 



* Figures in parenthesis indicate references listed following the text. 
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its optimal value before a significant change, A, in the objective func- 
tion occurs. The variation of the may be independent of the other 

or several parameters may be changed simultaneously in a co-ordinated fashion 
away from the optimum. The latter variation might allow for significantly 
large departures from the optimum for a specified A when the function 
possesses tne cnaractemstLCs or an X— cmensior4al ridge , 

In this respect, it is intended that this thesis may be used to evalu- 
ate and analyze the optimum of any general function of N independent 
variables in such a way that a complete sensitivity at the optimum is 
clearly presented in a useful output format. It is shown to the user 
which of the specified variables may be changed and the magnitudes of those 
changes before specified degration in the objective function will be 
exceeded. 
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II. PROBLEM DEFINITION 

An analysis of the behavior at the optimum of an N-variable function 

is possible where the second derivatives are known. Suppose that Y is 

a real valued function of N variables with continuous first and second 

partials and possesses a relative minimum at X , then the first deriva- 

o 

tive will vanish and by Taylor series expansion: 



Y-Y =AY=-;rAX 
opt 2 o 



>x 



AX^ + Higher Order Terms 



( 1 ) 



where 






is the Hessian matrix (H) of Y at the optimum. 



From the Taylor series approximation (1) we find that the gradient 



A Y(X) 




VY(X) 



and solving for X^ yields: 



( 2 ) 



X = X 
o 




VY(X) 



( 3 ) 



so that if 




were known, the step to the optimum would be 



given by (3) . Some optimization algorithms collect information which 



generate an approximation to the inverse Hessian matrix during the search 



for the optimum. These algorithms provide H ^ so that the problem of 



obtaining H, in the use of such an algorithm is reduced to the relatively 

■32 ‘ -1 

simple additional task of inverting the N x N matrix — ^ 

dx' 

Historically the method is similar to Newton’s method v/hich minimizes a 
function Y(x), x on N-vector, by generating a sequence of points (compare 
with equation (3)) 



^(k+l) ^ j^(k) _ ^(k) ^(k) 



where OL ' is an appropriately chosen scaler constant. Fletcher and 

Powell’s version is in fact a quasi-Newton method which uses an initial 

(k) - 1 

estimate and computational history to generate an estimate to [H^ 
rather than performing the computational work of evaluating and inverting 
the matrix at each step. 
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III. DAVIDON'S METHOD 

Davidon introduced a variable metric method which was the first 
optimization technique to use past information estimating the inverse 
Hessian matrix. Fletcher and Powell have improved upon the method by 
simplifying the iteration scheme and modifying the criterion for con- 
vergence. Their method, which numerically determines the minimum of 
a function of several variables, combines the best features of the con- 
ventional gradient method and Newton’s method, namely the sureness of 
convergence of the former and the quadratic terminal convergence of 
the latter. An excellent exposition of the method, including conver- 
gence proofs, has been given by Fletcher and Powell. For the purpose 
of this thesis however, it will suffice here to state the algorithm and 
point out the usefulness of its main features as was done by Kelley [3]. 

Let us denote the free parameters as and the function to be 

minimized as Y(x). It is assumed that Y has continuous first and 

o 

second partial derivatives with respect to X. Any starting point X 

o 

is chosen according to some criteria*. At the starting point X , the 
gradient vector Y as well as Y itself, is evaluated. A change is 

X 

then made in X^ according to: 

AX = -a Y (4) 



* The freedom in the choice of the starting point depends on Y . If Y 
is well behaved, this choice is free. In other cases, the starting 
point must be close to the optimum to assure convergence. 



where H is a positive definite, synunetric matrix, defining a metric 
in the X-hyperspace* Its initial selection is otherwise arbitrary. For 
general purposes the unit matrix may be used, but if the parameters differ 
by orders of magnitude it is convenient to either scale them or to estimate 
the accuracy with which they are to be determined, a > 0 is a step size 
parameter. 

In Davidon*s method, the one-dimensional minimum of F vs. a is 

O -1 

obtained along the vector originating in X in the direction 11 Y . 

X 

At the new X, the gradient vector Y is again evaluated. The 

X 

matrix is updated according to 

= H '^ + A + B (5) 



where 

*T 

AX AX 

A = 

AX AY 

X 



and generates the inverse of H in N steps for the quadratic function 
and where 



„ - H Ay ay 

B = XX 



H 






Ay H Ay 

X X 



is intended to cancel the initial assumption for H t ^ • The procedure 

* o 

is begun again with the new values of X, Y^, H 
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It is shown in [1] that H remains positive definite and that, as 

X approaches the minimizing point, approaches Y ^ evaluated at 

xx 

the minimum. For quadratic Y, in N-dimensional space, the minimum is 
obtained in, at most, N steps (within round-off error): the method is 

quadratically convergent. 

For more general functions having the smoothness properties assumed, 
a Taylor expansion through quadratic terms provides a good representation 
of the function in some neighborhood of the minimum. With converged, 

the minimum of Y vs. a then will occur for a = 1, 

AX = -a Y 

X 

will approach the value given by Newton’s method, namely 
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IV. SENSITIVITY OF THE PERFORMANCE FUNCTION AT THE OPTIMUM 

Assuming that the optimum and H at the optimum have been deter- 
mined and that H ^ has been inverted, then all the necessary information 
has been collected for a detailed analysis of the sensitivity of the per- 
formance function to changes in the performance parameters. 

Let the criterion function be 



=^2 V 



and assume for simplicity that the origin, i.e. 



X 



1 



= 



X = 0 

n 



is an analytic optimum of Y. Around the optimum 



Y 



1 

2 



x"^HX 



( 6 ) 



where 



H = 



3 2, 



9x . 9x . 
1 3 



= a. 



ij 



is the Hessian matrix of Y at 0. 



Independent Variation (one at a time) 

If one parameter at a time is allowed to depart from the optimum, 
the function Y is 



where a^^ is the corresponding diagonal term of H, For an assigned 
change A in Y we find 



2AY ' , ^ 

IX 

where ± Ax^ is the allowable change in for the previously assigned 

acceptable variation in Y, 

Simultaneous Variation 

If one allows several parameters to be changed in a co-ordinated fashion 
away from the optimum, departures from 0 for a specified Ay may be 
significantly larger than those shown in the one-parameter-at-a-time case 
To illustrate this, consider the case in a 2-dimensional parametric space 
where the function y presents the characteristics of a ridge as illustrated 
in Figure 1. 

The Y =1A line intersects the x^ and x^ axes at distances which 
are far less than the values of x^ and x^ at the end of the Y =1A 
contour. (Five times less in the case of Figure 1). 
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A Skew Ridge Illustrating Advantage of Simultaneous Variation 



Figure 1 
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Eigenvalues and Eigenvectors of H 

To be able to analyze the characteristics of the Y surface in the 
neighborhood of 0, it is convenient to make use of the eigenvalues and 
eigenvectors of H at that point. 

It is assumed here that the final use of this analysis will be in 
computer programs, and that matrices eigenvalue and eigenvector search 
subroutines are either available or easily implementable. To that respect, 
we note that H is a symmetric, positive definite matrix, and that know- 
ledge of this fact simplifies the implementation of such subroutines. 

Let and be respectively the eigenvalues and eigenvectors 

of H [ 5] , i.e. 



det (H - XI) = 0 

" H T. = X. V. 

1 11 

Since H is real and symmetric the eigenvalues are real and the 
eigenvectors are orthogonal and may be expressed in orthonormal form: 

1 (normal) 

0 i ^ j (orthogonal) 

If the N-parameters are allowed to be changed in the direction of the 
eigenvector, i.e. M V^, then the degradation in the perfor- 

mance function AY is: 




-T - 
V. • V. = 
1 3 
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( 8 ) 



AY 



1 - T 

^ k . V . H k . V . 

2 1 X 11 



but 



H V. = X. V. 
1 11 



and 



— T — 
V. • V. 
1 1 



1 



so that 



AY 






X. 

1 



and for an allowable variation A in Y 



(9) 



( 10 ) 

represents how far in the direction of the normalized i'’“ eigenvector 

one may travel on the response surface before degrading the performance 

function by AY. It is observed that k is a maximum for X . , 

min’ 




so that the direction of least 
sensitivity to changes in the 
parameters^ is given by 

the eigenvector with X = X . . 

Q 



In other words the length of the AX vector 
a given ‘ AY in the performance function. 




is maximized for 



Consider a hypothetical 2-variable optimization problem where the eigen- 
values of H at 0 assume the values of 1 and 10, i.e. 
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= 1 ; = 10 

and the associated orthonormal eigenvectors are: 



= (1,0) and = (0,1) 



respectively. The contour of the objective function at the optimum will 
then assume the shape of Figure 2. 

Observe that for a given Ay, the distance from the optimum in the 
direction of is considerably greater than in the direction of V^. 

In fact, it follows directly from (10) that the distance is/ ^ greater 



1 / 

It is evident that the relative length of any existing ridges in the 
response surface will be determined by the square root of the ratio of 
any two eigenvalues and that the magnitude of the sensitivity for a 
given A will be detennined by the square root of the eigenvalues. 

Let two distinct solutions of optimization problems assume the values 
Solution 1 (as before) Solution 2 



^11 ^ ^12 



^21 ^22 ^ 



= (1,0) = (0,1) 



= (1,0) = (0,1) 



The response surfaces will then assume the shapes of Figure 3a and Figure 
3b respectively. 

Notice that the ratio of the eigenvalues has remained constant and 
that the shape of the response surface is unchanged. The magnitude of 
the eigenvalues was decreased by a factor of 10 and if drawn to scale the 
response surface would be expanded in all directions by a factor of /lO . 
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Effect of Eigenvalue Ratio on Response Surface 



Figure 2 
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Figure 3a 



X2 




Effect of the Magnitude of Eigenvalues on the Shape of the 



Response Surface 



Figure 3 
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Skew Ridges 



The previous examples have pointed out that when the matrix H has 
one eigenvalue much smaller than all other ones a ridge will exist in the 
response surface. A skew ridge exists when one eigenvalue is much smaller 
than all other ones, and the corresponding eigenvector is not parallel to 
one of the axis in X-space. 



If a ridge is parallel to one 
axis in X-space , it can always 
be removed by a change in scale 
along that axis. It is there- 
fore representative of the way 
scales are chosen rather than 
a characteristic of the func- 
tion to be optimized. 



The situation is completely different if the ridge is at an angle to 
the axis because no change in scaling can remove the ridge. Such a ridge 
reflects an interaction between parameters in the way in which they affect 
the criterion function. 

Only when the eigenvalues and eigenvectors of the Hessian matrix are 

known, can such a ridge be discovered; and only then can the characteristics 

of the optimum be determined. complete search of eigenvalues and eigen- 

vectors of H and the derivation of the resulting sensitivity information 
is contained in the computer program described in the following section. 
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V. COMl'UTER PROGRAM 

The program described in this thesis makes use of a standard, 

Powel 1 -Fletcher type, parameter-optimization computer software 
package (6). It contains a multi-dimensional optimization algorithm 
similar to Davidon's variable metric method as was previously described. 

The matrix inversion and eigenvalue analysis subroutines were obtained 
from the scientific subroutine library of the IBM 360 model 91. A listing 
of these subroutines is presented in the Appendix. A simplified flow 
diagram illustrating the coordinated use of these programs is shown on 
the following page. 

Necessary Input for Sensitivity Output 

The input is the expression relating the objective function Y to 
the N independent variables x^. The expression may take the form of 
a single equation, e.g. 



Y = f(x^) 

or may comprise any number of subroutines as long as the objective function 
Y and the N independent variables are defined. The allowable departure A 
from the optimum is also a required input and must be specified by the user. 
It may be expressed as a percentage change in Y, e.g. 



Y = f(x) 






I 

I 

V 

POWELL - FLETCHER 
OPTIMIZATION ALGORITHM 



' ^OPT’ ^OPT’ ^ 



-1 



MATRIX INVERSION 



T 



H 



EIGENVALUES, 

EIGENVECTORS 



I H, A V 
V 1 ^ 



SENSITIVITY 

EQUATIONS 



AX. 

1 



(independent) 






(AX)^ (along eigenvector) 



Flow Chart of Computer Programs 
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Ay = 1% Y 



opt 



or Y = 99% Y 



opt 



or as a fixed quantity, e.g. 



AY = 10. 



Output Format 

Given the required input ; 



Y = f(x^) 

Y “ A (a constant) 



the sensitivity information is presented in the following table for con- 
venience in analysis. 



TABLE I 



Y ^ 
opt 


AY 


X . 

opt 


Ax fo- 

Independent 

Variation 


r Specified A 
Simultaneo 
^min 


y 

us Variati< 

X • • • • 


|[ 

1 


— 




^1 


±Ax^ 




• • • 


Axi 






^2 

1 


±Ax2 

1 


Ax^ 

1 


1 


Ax^ 

1 






1 


1 


i 


1 

i 


i 

% 
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VI NUMERICAL EXAMPLE: ROSENBROCK FUNCTION 

In order to illustrate the usefulness and accuracy of sensitivity 
analysis a standard test function is offered as an example. The function 



Y(x^, X2> = lOOCx^ - + (1 - 



was proposed by Rosenbrock and is interesting in that it possesses a 
steep-sided ridge following the curve = x^ as shown in Figure 4. 

The exact solution of the problem is offered along with the computed 
values so that a comparison can be made and the effectiveness of the 
techniques employed may be evaluated. 

Problem: minimize Y = 100 (x^ (1.0 - 

Exact Solution: 



Y = 0.0 
opt 





802 


-400 




0.5 


1.0 


H = 






II 








-400 


200 




1.0 


2.005 



X2 = 1.0 



where all numbers given are exact. 

Fletcher and Powell’s version of Davidon’s variable metric technique 
found an optimum solution: 

Y = 10~® - 0.0 
opt 

x^ = 1.00007 - 1.0 



X, 



= 1.00017 - 1.0 









I 



X2 




Rosenbrock's Curved Valley Function 



Figure 4 
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and the approximated inverse Hessian matrix 



H 



-1 



0.5055 

1.0051 



1.0051 

2.0035 



which represents an error of less than 1%. 

The matrix inversion subroutine with the approximated H as its 
input resulted in the Hessian matrix: 



H = 



828.3 



-415.5 



-415.5 

209.0 



in error of about 4%. This error is more than acceptable for the purposes 
of sensitivity analysis. 



TABLE II 



Sensitivity Data of Rosenbrock function: 









AX 


for Specified Ay 


Y ^ 


AY 


X. 

opt 


Indep endent 


Simultaneous Variation 


Opt 




Variation* 


X . = .399 

min 


X 2 = 1037. 


10'® 


0.10 


1.00007 


±.0155 


Ax = +.319 
^ / 


Ax^ = +.012 






1.00017 


±.0309 


Ax2 = +.633 


Ax2 = -.006 



Analysis : 

If the response surface were to be constructed from the sensitivity 
data shown with no knowledge of the function under consideration it would 
resemble Figure 5. 
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Predicted Shape of Rosenbrock Function at the Optimum 



Figure 5 
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A comparison with Figure 4 shows that the sensitivity analysis at 
the optimum gives the desired results. The above figure shows that the 
function is extemely sensitive to changes along the ^ixis = ±.015 

for Y = 0.1) as well as to changes along the axis (Ax^ = ±.031). In 
contrast, it is shown that a skew ridge exists along the direction defined 
by Ax^ = .318, Ax^ = .633 and that if varied simultaneously x^ may 
be changed by 30% and X 2 by 60% before Y changes by 0.1. 
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VII. EQUALITY CONSTRAINTS 

In the optimization of multi-variable problems it is often the case 
that only certain combinations of the variables are either meaningful 
or acceptable. The imposed restriction usually assumes the form of an 
equality (or inequality) constraint: 

G(x ) =0 
i 

The analytical solution of such a problem can be found by the method 
of Lagrangian Multipliers [7], which seeks values of the parameters for 
which the modified objective function 

F* = F + AG (11) 



is stationary, i.e. 



F* = f +XG = 0 

XXX 



( 12 ) 



Solving for X yields 



A 




(13) 



A is meaningful in that it represents the cost of the constraint G. If 
G were relaxed by 1 unit then F could vary by X. 



In the general case however, numerical search methods are unlikely 
to locate all types of stationary points for the modified function using 
Lagrangian Multipliers. A more feasible approach is that of penalty 
factors. [8] 

Penalty Factors 

In the use of penalty factors a modified function which incorporates 
the constraint is defined as 



where P is a large positive-valued constant (for minimization) . If P 
is chosen large enough then G is forced close to zero in the search for 
the optimum. At the optimal solution G is equal to some small quantity 
G. If e is not within an acceptable distance from G then P is in- 
creased until an optimum is found which satisfies G within e. The 
cost of the constraint can be found in a manner identical to the method 
of Lagrangian Multipliers. 

At the optimum of the modified function of (14) 



F* = F + PG^ 



(14) 



F * = F + 2Pe G = 0 

XX X 



(15) 



SO that 



2PC = -^=X 



( 16 ) 



X 



is the cost of the constaint G 
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Although the idea of a penalty function seems to have been conceived 
some years ago (Courant [1943]), there has been very little computational 
experience with regard to its applications. By examining equation (16) 
it is observed that as e approaches zero, P becomes infinitely large. 
Large values of P however effectively produce a sharp ridge in the con- 
tours of F*, and most search techniques are troubled by the existence 
of such a ridge. The choice of P is therefore a compromise between 
large values for small violations in G, and smaller values to eliminate 
troublesome ridges in the modified objective function. 

In using Davidon’s optimization technique however, it was discovered 
that even if a ridge presents no problem in finding the optimal solution, 

P may not be chosen arbitrarily large. For this situation, £ becomes 
so small that changes in the parameters of order e produce corresponding 
changes in F which are less than the criterion for convergence in the 
search for the optimum. If £ is to be meaningful it must be large enough 
to possess a unique solution. In other words changes in £ must be large 
enough to affect the terminal convergence of the optimization search. It 
is therefore necessary to possess some insight into the problem before 
the penalty factor method can be employed. 

We may note that for any given equality-constrained optimization 
problem,. X will possess a unique value. Analytically X is found to 
be 2P£ as £ approaches zero and P approaches infinity. Let us 

denote e as the maximum allowable violation in the equality constraint 
max 

and £ as the minimum £ which will produce a unique optimum within 

min 

the limits of the convergence criterion. £ must now satisfy 

£ . < £ < £ 
min max 

for a meaningful solution. 
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Values for X and P are estimated within one or two orders of 



magnitude such that 



2P f v£ — X. .. 

(est) (est) 

If the resulting optimum possesses an allowable e then the solution is 
found with the cost of the constraint 

X = 2Pe (17) 

and the constraint violation: 



e = G - G 
opt 



(18) 



If c > e or e < e . then P must be increased or decreased re- 
max min » 

spectively until an allowable e is found. 

Example Problem 

To illustrate the use of the penalty factor method for equality con- 
strained optimization problems the following two examples are offered for 
comparison and analysis. 



Let (x^ - 3)^ + (x^ - 3)^ 

and ^2^ = x^ + X 2 ~ 4 = 0 

G (x^, x^) = x^ x^ - 4 = 0 
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Figure 6b 



Interpretation of Cost of the Constraint 



-F 



X 



G 

X 



Figure 6 
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The response surfaces are shown in Figure 6, where the equality 
constrained function has a minimum value at the same point for both 
problems, i.e. = X 2 = 2.0. Since the unrestricted function is the 

same and the solution lies at the same point, then the gradient F is 

X 

the same for both problems. The constraint and its gradient are 

different however and consequently the cost of the constraint A will 
have two distinct values. 

Analytically it can be shown that -F /G = A = 2 for the linear 

X X 

constraint x + x^ = 4 and that -F /G = A = 1 for the hyperbolic 
constraint x^ x^ = 4. The cost of the constraint has been reduced in 
the second case because F is less sensitive to changes in G^ at the 
optimum. 



Computer Results: 



For 


F* = F + P 


and 


p* 


and 


P = 100 




p 




F* = 1.980 




F* 




x^ = 2.005 




^1 




x^ = 2.005 




^2 



£ = G* - G = .010 £ 

X = 1.988 X 



= F + P 

= 100 
= 1.994 
= 2.0009 
= 2.0020 
= .0057 

= 1.141 



In the second case, a minimum was found which was very close to the 
constraint, e.g. £ = .0057. As a result X = 1.141, a difference of 
.141 from the analytical value. This error may be explained by examining 
the magnitude of £. An £ of .0057 is very small for the problem under 
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consideration. Changes in the variables of order £ will not affect the 



optimal solution, so that £ is actually less than previously 

defined. 

A second optimization was performed with P reduced to 50 in order 
to increase £. (Remember that it has been shown that 2Pe is constant). 



The results were: 




F* 


1.989 


II 


2.0065 


X£ 


1.999 


e 


.011 


X 


1.07 



It is observed that decreasing P resulted in a more meaningful value 
for £ and consequently a closer approximation to the const constraint, X 
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VIII. SENSITIVITY ANALYSIS OF A NUCLEAR ROCKET ENGINE 

As was previously mentioned, the primary purpose of a sensitivity 
analysis is to aid the design engineer in constructing a system which 
will operate in an optimal fashion, despite variations in the controlling 
parameters. To that effect, an illustrative example is presented whereby 
the engine design parameters of a nuclear rocket are optimized in order 
to achieve a maximum payload at a specified hyperbolic velocity. 

A set of mathematical models of the elements of nuclear rocket 
engines, suitable for detailed systems analysis, has been developed [9] 
which constitutes the basis for a digital computer program called 
NUROC/SAC (Nuclear Rock et _^stems Analysis ^ode) . The Code has been used 
to describe a number of existing engines and the results obtained were 
found to be accurate [10] • 

ESCAPE [11] is another computer code which calculates geocentric 
(or planetocentric) tangential thrust escape trajectories and which may 
be used in conjunction with NUROC/SAC. 

Input Parameters 

From a design viewpoint the most important input parameters to 
NUROC/SAC are: 



Q 



thermal power of the reactor, watts 



D 



diameter of the reactor core, meters 



L 



length of core, meters 



L/D 



ratio of core length to diameter 



T 



maximum allowable core material temperature, 



cmax 




I 



»c 


nozzle chamber stagnation pressure, n/cm^ 


NAR 


nozzle area ratio: A . /A , 

exit throat 



The output format of NUROC/SAC consists of a sunmiary of operating variables 
including the input design parameters, defining entirely the characteristics 
of a specific nuclear rocket engine. The most important of these and the 
ones which will be used as inputs into ESCAPE are: 



* 

m 


total engine propellant mass flow 


F 


total engine thrust 


m 

E 


total engine mass 


Additional inputs 


to ESCAPE which must be specified are: 


m 

0 


initial mass in earth orbit 


H 

0 


initial orbital altitude 




final hyperbolic excess velocity specified 




tankage to propellant mass ratio 



The output of ESCAPE may assume a variety of formats but the payload 
delivered at the specified hyperbolic velocity is the payoff in the opti- 
mization problems which follow. The payload is defined as the initial mass 
minus the engine mass, the tankage mass, and the propellant mass required 

to reach V, 
hf 

A previous study [12] has indicated that within the range described by 
technological constraints the payload perfomance of the nuclear rocket 
will increase monotonically with the maximum allowable core material 
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temperature and the nozzle area ratio, and will decrease mono tonically 
with the chamber stagnation pressure. Thus choice of these parameters 
is limited to technologically realizable values. Typical values presently 
in use are: 

T = 3000 

cmax 

p^ = 300 n/cm^ 

NAR = 100. 

The power, power density, diameter, and length which parameterize the 
core geometry may be varied in a constrained but optimal fashion to describe 
an engine which will inject the maximum payload into a specified inter- 
planetary trajectory for a given initial mass in earth orbit. The problems 
which follow will consider an initial mass of 100,000 kilograms in an earth 
orbit of 500 kilometers and will optimize the core geometry to inject the 
maximum payload into a trajectory described by a hyperbolic excess velocity 
of 10,000 m/sec. 

Optimization of Core Geometry 

It is interesting to note that if any three of the four engine 
parameters, power (Q) , power density (^) , diameter (D) and length (L) , 
are specified then the fourth is automatically determined. Since the 
length and diameter describe the volume, the power density can be ex- 
pressed as a function of Q, L, D, i.e., 

2 

ip = AQ/ttD L 

Thus the optimization of core geometry is reduced to a 3-dimensional 
problem. 
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Maximize Payload: Free Parameters Q, L/D 

The following example maximizes the payload delivered by a nuclear 
rocket with free parameters: power (Q) , power density (^) y and ratio 

of length to diameter (L/D). 

Starting Point; 

Q = = 2000. MWt 

4^ = x^ = 6.0 X 10^ W/m^ 

L/D = ^3 = ^-0 



TABLE III 



Sensitivity Data at the Optimum 



Y ^ 
opt 


AY = 1% 


X ] 

opt 


j Ax f 

[nde pendent 
Variation 


:or Specifie 

Simultaneo 
X = 39.5 


:d Ay 

us Variatior 
X = 1105. 


1 

X = 6150 


1,630 Kg. 


316. Kg 


Q=1832 MWt 


AQ=±333. 


68. 


-28. 


28. 






4^=8.84x10® 


A4^=±3.03 Q 


3.88 


-.062 g 


-.017 






W/m® 


(xlO^) 


(xlO^ 


1 (xlO^) 


o 

r— < 

X 




L 


/D=4.80 


AL/D=±.57 


.387 


.677 


.199 



Analysis : 

When both ip and L/D are free parameters the optimized values 
become so high that the cost of the uranium necessary to make such a 
reactor critical would become prohibitively expensive. Also, a length 
to diameter ratio of approximately 5 and a power density of 9x10^ W/m^ 
would describe a highly inefficient reactor due to excessive neutron 
leakage through the core ends. Other important design factors (such as 
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shielding), which are not now contained in NUROC/SAC, indicate that the length 
should be only slightly larger than the diameter. It should also be noted 
in the sensitivity data that ip is quite flexible (A^ = ±3x10^) and 
that smaller values may be chosen with little resulting loss in payload. 

This problem should serve as an example that one cannot randomly 
optimize variables or undertake a sensitivity analysis without first 
acquiring a knowledge of the system under consideration. With these 
facts in mind, a second optimization problem is solved in two-dimensions 
with L/D fixed at 1.5. 

Maximize Payload: Free Parameters Q, 

Starting Point: 

Q = = 2000 MWt 

ijj - - 6.0x10^ W/m^ 

"TABLE IV 

Sensitivity Data at the Optimum 









Ax for Specified Ay 


Y 


= 1% X 


Independent 


Simultaneous Variation 


opt 




opt- 


Variation 


1 X = 5151 


X = 23630 


30,938Kg. 


309. 


Q=1823 MWt 


AQ = ±339. 


340. 


9.6 






4^=4.97x10^ 


Aijj “ ±.160 


-.02 


.159 






W/m^ 




9 


9 








(xlO^ 


(xlO^) 


(xlO ) 



Analysis : 



The restriction on L/D (fixed at 1.5) results in a loss in optimal 
payload of only 2 percent. It is interesting to note that the optimal 
power and allowable deviation are almost identical to the three-variable 
problem. The power density however has a considerably lower optimal 
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value and is much smaller when L/D is fixed. 

The sensitivity data reveals that any power between 1300 and 2150 MWt 
may be used if is held constant with only a small (1%) resulting loss 

in payload. In contrast the power density must remain close (v/ithin 4%) 
to 5.0x10^ W/M^ . The fact that the eigenvectors are located close to the 
Q and axes implies that little more can be gained by varying Q and 

simultaneously. 

Maximize Payload: Free Parameters Q, L/D 

With the discovery that and L/D cannot both be free and having 

fixed L/D at 1.5 it would now be advantageous to fix and optimize 
on Q and L/D. A value of 3.0x10 W/M is chosen for the power density 
based on accepted values for existing nuclear rocket engines. 

Starting Point: 

Q = f 2000 MWt 
L/D = ^2 = 2.0 



TABLE V 



Sensitivity Data a t the Optimum 



%t 


AY = 1% 


X 

opt 


Ax for 
Independent 
Variation 


Specified Ay 
Simultaneous 
X = 935. 


Variation 
X = 3867. 


30,675 Kg. 


306. 


Q = 1542 MWt 


AQ = ±397 


AQ = 114 


AQ = 390 






L/D = 1.67 


AL/D = ±.78 


AL/D =0.79 


AL/D = -.0561 

1 

1 



Analysis : 

When the power density is fixed both the power and the ratio of length 
to diameter have large acceptable variations. The power may vary from 1150 
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to 1950 with L/D and ip fixed, and L/D may vary from 1.9 to 2.45 with 
Q and ip fixed. This allows for considerable flexibility in designing 
an optimal engine as long as the power density remains relatively constant. 

The least sensitive direction of change is associated with the mini- 
mum eigenvalue. The associated eigenvector indicates that L/D may still 
vary from .9 to 2.45 and Q need not remain fixed, but may vary by ±114 MWt 
as long as the direction of change from the optimal value coincides with 
changes in L/D. 

Verification of Results 

In order to verify that the predicted sensitivity of the optimized 
variables is accurate, the allowable deviations were substituted into 
NUROC/SAC and ESCAPE and the resulting payload was computed. A comparison 
was made to check if the payload remained within the predicted 1% of the 
maximum. The optimal values were first rounded-off to Q = 1550 ±400 MWt; 
L/D = 1.67 ±.75. 

TABLE VI 

Comparison ; 



Maximum payload minus 


1% = 


30,370 Kg. 




Independent Variation 




L/D 


Payload <cKg) 


I. AL/D, Q fixed 




.90 


30,271 






2.45 


30,520 


II. AQ, L/D fixed 












1150. MWt 


30,295 






1950. MWt 


30,420 


Simultaneous Variation 




L/D 






1435 


.90 


30,135 




1665 


2.45 


30,445 
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The table of comparisons indicates that if Q and L/D are increased 
either independently or simultaneously the loss in payload is even less 
than the predicted 1%. If decreased the deviations become slightly larger 
than 1% but are still highly accurate. This may not always be the case. 

It must be realized that the sensitivity analysis is accurate only 
at the optimuni. For functions which are very flat (i.e. insensitive to 
change) the predicted variations may be quite large, as indeed they are 
in the example given. When the variations are substituted, the function 
is no longer close to the optimum and results may vary considerably from 
those predicted. For this reason it is prudent to verify results especially 
at points far from the optimum. 

Another reason for verifying results is that there exist small in- 
herent errors in the optimization search, matrix inversion and eigen 
analysis subroutines. Care should also be taken in defining the conver- 
gence criterion because if. the computer is forced to make repeated searches 
near the optimum the values for the inverse Hessian matrix will be destroyed. 
Any resulting sensitivity analysis will be meaningless. 

Application of Penalty Factor 

In the previous two-dimensional optimization problem an optimal power 
of 1550 MWt and length to diameter ratio of 1.67 was established for a 
fixed power density of 3.0x10^ M/W^. As was previously mentioned for 
equality constrained problems it is of interest to know the cost of the 
restriction on power density X. The following example employed the 
penalty factor method to optimize a three-variable problem with the constraint 

G(x^) = - 3.0x10® = 0 
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Starting Point; 



Q = = 2000 MWt 

L/D = ^2 = 6.0 

= 4.0x10^ W/M^ 



Penalty Factor = P = 10,000 



Results : 



Payload = 30,670 Kg. 

Q = 1547 MWt. 

L/D = 1.7 

llJ = 3.020x10® W/M® 

£ = G* - G = .020 

X = 2P£ = 406. 



The cost of keeping the power density constant is 406 Kg. of payload 

9 

per (10 W/m^) . 



-42 
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IX. CONCLUSIONS AND RECOMMENDATIONS 

A sensitivity analysis at the optimum has been performed through 
the use of existing computer codes. By simply defining an objective 
function and any number of independent variables an optimized solution 
is found. In addition, it has been shown that the sensitivity of each 
of the variables at the optimum may be given in a useful format. With 
all others held constant, the range that each variable may assume 
before a specified degradation in performance occurs is given. Also, 
the length and direction of the axis of all ridges in the objective 
function are listed in order of decreasing magnitude. With this 
information at hand a complete knowledge of the sensitivity of all the 
input parameters is available to the user and decisions regarding 
optimal parameter choice can be made, taking into account the 
flexibility given by the sensitivity knowledge. 

The methods employed were shown to be highly accurate when dealing 

with purely mathematical problems. The sensitivities of Rosenbrock's 

parabolic valley function were found with little difficulty. Engineering 

applications, on the other hand, require a great deal of insight into 

the problem before a sensitivity analysis may be started. Scaling the 

variables is important so that the sensitivity information is meaningful. 

The nuclear rocket example studied the engine performance related to 

power, power density and ratio of length to diameter of the nuclear core. 

A scaling problem existed because the power in watts and power density 

9 

in watts/cubic meter were of the order of magnitude of 10 , while the 
length/diameter was approximately unity. Usually such a problem can be 



avoided by normalizing each variable by dividing it by the order of 
magnitude. This would be fine if the acceptable deviations occurred 
within the same number of significant figures. However, if one variable 
is much less sensitive to change than all the others, then the eigen- 
vector associated with the minimum eigenvalue will be dominated by that 
component and very little knowledge can be gained about the other 
variables in that direction. It would be convenient to modify the 
computer program so that the sensitivities would be normalized instead 
of the variables. Of course, the user could always eliminate the 
problem by normalizing the variables the first time, and after observing 
the resulting deviations, normalize the sensitivities and optimize 
again. The optimized solution will be the same but the sensitivity 
data will be more accurate and meaningful. 

Another problem experienced when working with engineering problems 
was in defining a convergence criterion. For the Rosenbrock function a 
change in the variables of 10"^ produced significant changes in the 
objective function. When optimizing the nuclear rocker, however, the 
maximum payload was essentially determined when the normalized variables 
were accurate in the second decimal place. With the convergence criterion 
set at 10“^, the program made over a hundred more iterations before 
stopping and found an optimum which was only about one kilogram more in 
payload. When the optimization search stays close to the optimum for 
many iterations, the inverse Hessian matrix is destroyed and any 
resulting sensitivity analysis has no meaning. A modification in the 
program is required so that when the optimum has essentially been 
determined, the inverse Hessian matrix can be stored and any subsequent 
refinement of the optimum will not affect the sensitivity data. 
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Equality constrained optimization has been discussed but not in 
detail. The penalty factor method used by Kelley was implemented 
to determine an approximation to what is usually referred to as the 
cost of the constraint. 

It would also be of interest to know how far in the plane of the 
constraint and normal to it, one may travel on the response surface 
before reaching a specified change in the optimum. This thesis has not 
considered such a sensitivity analysis. The techniques involved are 
similar, but more attention is needed in this area. 

Inequality constrained optimization problems are generally no more 
difficult to solve than equality constrained ones. For well behaved 
functions, the solution either does not violate the constraint or lies 
on the boundary and may be treated as an equality constraint. In 
either case, the techniques developed in this thesis can be applied. 
Care should be taken to try, a variety of starting points so that if the 
same solution is reached, the function can be assumed to be well 
behaved . 
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appendix a 

This appendix describes the implementation of the variable metric op- 
timization method and matrix inversion and eigenvalue analysis computer 
programs for the purpose of sensitivity analysis. The input variables 
and their definitions are shown along with a flow diagram illustrating 
how the programs were modified to run successively and produce the desired 
sensitivity data. A typical output format is also given. 

The matrix inversion and eigenvalue analysis subroutine listings con- 
tain brief descriptions of the methods employed and their references. For 
a complete description of the computer coded variable metric method consult 



reference 6. 



fj 






BLOCK DIAGRAM OF CONTROLLING PROGRAM 



CALL VARI-IET 



START 



(N, F, X, XXIN, XMAX, 
DELX, STEP, a, B, 
NEXP, della:-:, change, 

JPOINT, JFLOW, H) 



1 LOAD 

I N, STEP, a, 3, 
i NEXP, DELLAM 
i CHANGE, JPOINT 
X(J), X(J)MIN, 

. X(J)MAX, DELX(J) 



j 

i 


TEST 


FALSE 1 1 


N < 


10 ^ 


HALT 

f . . - 


X(J) 


> XMIN(J) i 


FALSE XMINIJ')=X('J') 


|X(J) 

1 


< XMAX(J) * 


XMAX(J)=X(J) 




j TRUE 


-1 



RETURN 



/\ 




PRINT OPTIMUM AND 

I 

CO-ORDINATES 
i FB, X(J) J=l, N 

I 

TEST CONSTRAINT 
G = 0 



TRUE 

! PRINT INVERSE 
HESSIAN PLATRIX 
H(I,J) 

INITIAL CALL , ' 

MODEL (N,F,X) j 

CHANGE STORAGE 
! MODE 

HH(I) = H(I,J) 



NEVAL = 1 
FB = F 
JFLOW = 0 



NGRD = 0 
NIMP = 0 



INITIAL PRINTOUT 
NIMP, NGRD, NEVAL, 
FB, X(J), J=l, N 



MATRIX 

INVERSION 



OPTIMIZATION LOOP 



( 






r 

\ 



CALL MINV 

j (HH, N, D, L, M) 
' RETURN 

PRINT HESSIAN 
MATRIX 
HH(I) 



CHANGE STORAGE 
MODE 

H(I,J) = HH(I) 

^ ^ 

STORE DIAGONALS j 

UNCOR(I) = H(I,I) 



i I 

i 

I 



EIGENVALUE 

ANALYSIS 



CALL EIGEN 

i 

•>( 

I (HH, R, N, MV) 

i ~ 

; RETURN 

t 

I GRANGE STORAGE MODE 
I H(I,J) = HH(I) 



f 

: STORE EIGENVALUES 
^ DECEIG(I) = H(I,I) 



j SENSITIVITY EQUATIONS 
i 

I DY = 




■ V. = k. V. 

. 1 XI 



! 



V 

READY FOR OUTPUT 




PRINT 

OPTIMAL Y VECTOR OF OPTIMAL X 

ASSIGNED VARIATION 'DY' 

INDEPENDENT VARIATION ±DX. 

1 

SIMULTANEOUS VARIATION X., V.. j 

I 

__i 



V 

EXIT 



Variable Metric Program 



Input variables and definitions: 



N 


= the number of independent variables 


STEP 


= initial step size for X components 


ALPHA 


= factor by which step size is increased (a > 1) 


BETA 


- factor by which step size is decreased (0 < 3 1) 


NEXP 


= limit on number of experiments along a vector line 


DELLAM 


= termination criterion for vector search expressed as a 

range fraction 


CHANGE 


= program termination criterion expressed as a range fraction 


JPOINT 


= control for type of numerical differentiation desired in 

gradient subroutine . 

1 = forward difference 

2 = backward difference 

3 = central difference 


X(J) 


= vector of initial values for independent variables 


XMIN(J) 


= lower bound for each independent variable 


XMAX(J) 


= bound for each independent variable 


DECX(J) 


= step size for each independent variable used in gradient 

subroutine for numerical differentiation. 



There are N + 2 data cards required as inputs for each program execution. 
The first is an identity card written in any format. The second card contains 
the first eight input variables listed above in the following format fields. 
(110, 3F10.4, no, 2F10.4, 110) 

Each of the N remaining data cards contains the initial value, range 
and step size of the independent variables according to the following format 



statement. 


FORMAT (4F15.5) 
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I 





I 



r 



1 



l! 



A typical output frot: the variable metric optimizauion pregram, follows. 
It should be noted that for the inverse Hessian matrix to be printed properly 
the array designated H must be dimensioned exactly > i.e. if N = 2 then 
DIMENSION H (2,2) • 



A-5 



J 



)B0 



2 



:ep 




0.1 0000 


.PHA 




1. 00000 


•TA 


- 


0.50000 


:xp 




5 


!LLAK 


- 


0.01000 


lANGE 


= 


0. 0C01C 


>CINT 




1 



C. C010C 

0.00 100 



XHIN 



XMAX 



- 1.20000 
X- 1.00000 — 



•5.00000 

•5.00000 



5. 00000 
5. COOOO 



PENDENT VARIABLES 



0. 12000E 01 C. 1C0CCE 01 
C. 51186? 00- -0.35763E-06 
C. 35741E 00 -C. 2244i<E 00 
0. 14259E 00 0. 89546E-01 
C.5S605E-07 -0.41723E-06 
0.62192E-06 C.11936E-C7 
C.20230E-06 -0. 31673E-08 



TIMIZATICN COMPLETE AFTER ■ ■ ^3 FUNCTION ' E V ALU AT 10 K’S 

6 9 53 0.39664E-13 -0.20233E-06 -0.325962-08 

THE INVERSE HESSIAN MATRIX IS 




C.5CC1 



HESSIAN MATRIX 

' W A 9590 “fe:. ^ 

■^‘■■'"^■-" 1.5983 — 

- 1.9983 
3.9973 



*ERRCR**» ^ EN-1 

PRCGBAMME WAS EXECUTING LINS 65 IN ROUTINE M/PROG WHEN TERMINATION OCC 



NIMF 


NGRB 


NEVAL 


FUNCTION INDE 
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0. 10833E-01 
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C.40146E-12 
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Matrix Inversion and Eigenvalue Analysis 

The example output of the variable metric computer code shov/n previously 
lists the optimum point and the inverse Hessian matrix at the optimum. 

Before the matrix inversion subroutine (MINV) may be implemented the upper 
triangular elements of the symmetric H matrix must be stored in a singly 
dimensioned array (HH) • This was necessary because of the input format 
utilized in the calling sequence of MINV and was accomplished by inserting 
a ’DO’ loop in the main program. 

Once the matrix has been inverted, the elements of the Hessian are 
printed and the diagonal elements stored for calculation of sensitivity 
information. Implementation of the eigenvalue analysis subroutine 
(EIGEN) merely requires defining the constant MV, i.e. 

MV = 0 eigenvalues only 

MV = 1 eigenvalues and eigenvectors 

A dimension statement in the main program for sufficient storage space of 
the vectors L, M, (utilized by EIGEN) and R (utilized by MINV) is also 
necessary. All of the required information to construct the table of 
sensitivity data is now available in the main program. Insertion of the 
proper equations and output statements completes the modification for 
sensitivity analysis. 

A listing of MINV and EIGEN and an example output format follows. 



I 
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SOLMML'T :n 1: .'IKJV 

F-Ul^rosi: 

INVliHT A ^;ATI^'IX 



USAGE 

CALL .*^1 NV (A, N,D,L/i-i) 



USSCKIEIION OF FA i< Al-i E 1 F F< S 

A - I X P U T X A T K I X , D E S 'I S C Y E D IX C C Xi I U I 
FE SUL I A XT IXVEPSE. 

X - CXCER 0? X ATP IX A 
D - KESULTAXT C ET E PM I X A X T 
L - XOSK VECTOR OF LEXG"f'E X 
iX - X'CBK VECTOR CF LENGTH X 



10 X AXj 



. i' i. 



'I 



REMARKS 

MATRIX 



MUST OE A GENERAL MATRIX 



SUBROUTINES 

NONE 



AND FUNCTION SUEPKOC-SAMS REQUIRED 



THE STANDARD GAUSS-JORDAN METECD IS USED. THE DE 
IS ALSO CALCULATED. A DETERMINANT OF ZERO IXDICA 
THE MATRIX IS SINGCLAR, 



i r.i 
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SU’JbCU'IlM- Ml .W (A, \ 

A { 1 ) , L ( 1 ) , M (1 ) 



15 A DOUIILF :'HF:CI.SI0N V 
C IM CCLWMM 1 SHC:iLE B£ 
STATEJl-NI WlilCH FOLIC'^S 



US ION C1‘ TI: IS S CL' TIN 
SEMCV::D Tlil DCL 



IS 



15 I- 



DCUBLF FFECI5ICN A , D , b' I G A , HC L D 

THE C EUSI ALSO BE R2ECVED EKCE COLDLE PR5CISICN 
APPEARING IN OTHER SCUTINES USEE IN CONJUNCTION 
RCUTIKE. 



*-■ '> / r' - 

^ I ii. A ^ 

/. T T H C 




THE DOUELC PRECISION VERSION CP THIS 3U5RCUTI 
CONTAIN EOUBL’E PRECISION FORTRAN FUNCTIONS, 

10 MUST 3E CHANGED TO CAES. 



NE NUST ALSO 
AES IN STATIEENT 



SEARCH ECS LARGEST ELEMENT 

D= 1 . C 
NK=-K 

CO 8C K=1,N 
NK=NK+K 
L(K)=K 
K (K) ='R 
KK = NN<-K 
3IGA=A (KK) 

DC 20 J=K,N 
IZ=N« (J-1) 

DC 20 I=K^N 
IJ=I2+I 

Ir{ ADS(3IGA)“ ABS(A{IJ))) 15,20,20 
EIGA= A (10) 

L(K)=I 
K (K) =J 
CONTINUE 

INTERCHANGE SCWS 



J=L (K) 

IP(J-K) 35,35,25 

XI=K-N 

DC 30 1=1, N 

KI=XI+N 

HOLC = -A (K3) 

JI=KI-K+J 

A (KI)=A (JI) 

A{JI) =HCIC 



INTERCHANGE COLUMNS 



I=M (K) 

lE(I-K) 45,45,38 
JP=N=* (1-1) 

DO 40 0=1, N 
JK=NK+J 

A-9 
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ll I J I* + J 
t!ClD--A 
A (Ji\) = A (JI) 

A(JI) =HCLC 

DIVICZ COLUMN 3Y MINUS IIVCI (VALUZ OP i^lVOl r I Z NT 10 
CONTAINED IN' EIGA) 

IF(3rGA) 48,46,48 
C= 0 - 0 
FETUHN 
EO 55 1=1, N 

50,55,50 

IK=NK+I 

A (IK) = A (IK)/ (-BIGA) 

CONTINUE 

HEDUCH MATEIX 

DC 65 1 = 1 , N 
IK=NK+I 
HCLD = A (IK) 

IJ=I-N 
CO c5 J=1,N 
10 = IJ +N 

If (I-K) 6C,c5,6C 

Ir(J-K) 62,65,62 
KJ=IJ-I+K 

A (10) =iiCLD=*A (KJ) +A (IJ) 

CONTINUE 



EIVICE soil BY PIVOT 

KO=K-N 

DO 7 5 0= 1 , K 

KO=KJ+N 

I?(J-K) 7C,75,7C 
A (KJ)=A (KJ) /BIGA 
CONTINUE 



PBOCOCT OF PIVOTS 
D=D“)‘EIGA 

RKFLACE PIVOT BY BECIPROCAL 

A (KK) = 1. 0/2IGA 
CONTINUE 

FINAL ROW AND COLUMN INTERCHANGE 

K = N 

K= (K-1) 

IE(K) 150,150, 105 
I=I (K) 

IF(I-K) 12C,120,108 
JC=N* (K-1 ) 

OR = S’* (I- 1) 

DO 110 J=1 ,N 
JK=JO+J 
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;iO L D - A ( J S ) 

11 = J i; + j 
\(JK)=-A (.11) 
i(JI) =IiCLD 
[=:•; (K) 

r (J-K) ICC, TOO, 125 
I = K-\' 

0 130 1=1, N 
I=KI+N 
iOLC=A (KI) 
ir=Ki-K<-j 

1 (KI)=-A (JI) 

(JI) =HC1E 

iC 10 ICC 
ETUBN 

:ne 
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SCBi^CUriM- 5IGZK 



PUKP0S3 

CC^'.r-U'IB EIGENVALUES ANC LIGZNVLCTCRS C- 
E AT H I X 



C V V ; ■ 



USAGE 

CALL EIGEN (A,R,N,MV) 



DESCRIPTICN C? FARAZETE5S 

A - ORIGINAL MATRIX (SYMMETRIC), DESTRCYEC 
RESULTANT EIGENVALUES ARE DEVELCFED IN 



IN CCMFjTAIIC: 
D I A GO N A i- 0 * 



MATRIX A IN I ESC ENDING CrUB. 

R - FESULTANT MATRIX OF EIGENVECTORS (STORED COLU M L N IS , 
IN SAME 3ECUSNCE AS EIGENVALUES) 

N - ORDER C? MATRICES A AND ?. 

MV- INFCT CODE 

0 COMPUTE EIGENVALIES AND EIGENVECTORS 

1 COMPUTE EIGENVAiUES ONLY (R NERD NCT 11. 
DIKENSICNED EUT MUSI STILL APPEAR IN CALLIN 
SEQUENCE) 

REMARKS 

ORIGINAL MATRIX A MUST 31 REAL SYMMETRIC (STORAS.'. MODE-V) 
MATRIX A CANNOT EE IN THE SAME ICCATICN AS MATRIX F 

SUBROUTINES AND FUNCTION SUEPBOG5AMS REQUIRED 
NONE 



METHOD 



DIAGCNAII2ATICN METHOD CEIGTNATED BY JACOBI AND ADA? 
BY VON NEUMANN FOR LARGE CCMPLTEHS AS FCUND IN 'MATH 
METHODS FOR DIGITAL COMPUTERS', EDITED EY A. RALSTON 
K.S. KILE, JOHN KILEY AND SONS, NEv; YORK, *1962, CHAP 
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SU 3 !' CU T I N E E .I' 2 N ( A , 'E , N , I'l V ) 
DIMENSION A(l),;ni) 



IE A DCU3LL' PHECISICN VE5SICN 0? 
C IN CCIUMN 1 SHCULC 31 SEI'.CVED 
STATEMENT isfilCJl FOLLOWS. 



T h' I S P 0 L T I N :• I : 
5 5CM THE DCEciLH 



DOUBLE 

1 



PRECISION 



A , K , A N 0 3 1 ! , A N R M X , -I H R , X , Y , S I M X , S I N N 2 , C C 3 , 
CCSX2, SINGS, FAKGI- 



THE C MUST ALSO 3?. 
APPEARING IN CTHEE 
SOUTINE . 



REMOVED ERCM DGUBL:: PEECISICN 3 
ROUTINES USED IN CONJUNCTION w I 







THE DOUELE PRECISION VERSION OF THIS SU3RCUTINE MUST 
CONTAIN DOUBLE PHECISICN FORTRAN FUNCTIONS. SCSI IN 
40, 68, 75, AND 78 MUST EE CHANGED TO DSCRT. A3S IN 
62 MUSI 3E CHANGED TC DABS. THE CCNSTANT IN STATEMTN 
3Z CHANGED TC UOD-12. 






5 SHOULD 



GENERATE IDENTITY MATRIX 

RANG E=1 . CE-6 
:IF(MV-1) 10,25,10 

IC = -N 

DO 20 J= 1 , N 
I0=IC+N 
DO 2G 1=1, N 
IJ=IC^1 
R(IJ)=C.C 
IF(I-J) 20,15,20 
B (IJ) = 1 .0 
CONTINUE 

COMPUTE INITIAL AN! FINAL NCRHS (ANORK AND ANORMX) 

ANORM=G.C 

DC 35 1=1, N 

DO 35 J=I,N 

IF (I-J) 30,35,30 

I A=I+ (J*J-J)/2 

ANORM = ANCHM +A (lA) (I A) 

CONTINUE 

I? (A NORM) 16 5,165,40 
AN03H=1 .4 14*SQRT (ANORK) 

ANRMX=ANOBM*RANGE/FLCAT (N) 



INITIALIZE INDICATCSS AND CCMEUTE THRESHOLD, THR 

IND=C 
THR=AN03M 
TKR=TrI3/FLCAT (N) 

L= 1 
M=L + 1 



COMPUTE SIN AND COS 



A-13 



in u. 



.'’.C= 'I) /2 

L0= (K-L-n/^' 

L!''.= I + ^C 

II-' ( .^G£ (A (LK) ) -THS) 13C,65,65 

TND= 1 

I.L=L+LC 

M ii = M -fr M C 

X = 0.•5^' (A (LL) - A (MX) ) 

Y = -A (LX)/ SCUr (A (IX) >)‘A (LX) +X^X) 

IF(X) 70,75,75 
Y=- Y 

SINX = Y/ SC r-T ( i. O'-" ( 1. 0+ ( SGRI ( 1. 0- y* Y) ) ) ) 

SIN X2=S IN X^i'S IN X 
CCSX= SCRI ( 1. 0-SINX2) 

CCSX2=C0SX’i-C0SX 
sixes -SIXX’»‘CCSX 

RC-rATii I AND X COLUMNS 

ILQ=X* (1-1) 

IXC=N>'‘ (X- 1) 

E-0 12 5 1= 1 , N 
IQ= (I>f‘I-I)/2 
Ir(I-L) 8C,115,80 
I?(I-X) 65,115,9C 
IX=I+iMC 
GO TO 95 
I« = X-HC 

IF(I-L) 1CC,1C5,1C5 
IL = ItLC 
GO 10 n C 
II=L+IQ 

X = A (IL) «COSX-A (IH) vsiNX 
A (IX) =A (II) *SINX + A (IK) *CCSX 
A (II) =X 

IF(XV-I) 12C, 125,120 

ILR=ILC+I 

lMR=lXQ-t-I 

X = R (ILR) *CCSX-E (IXR) *SINX 
S(IXR)=R(IlR)>i‘SIKX + R(IKR)*CCSX 
3 (I IB) =X 
CONTINUE 

X=2. 0>!<A (LX) ’i'SINCS 

Y=A (LL) =i‘COSX2 +A (XX) *SINX2-X 

X=A (LL) *SINX2+A (XX) *COSX2 + X 

A (LX) = (A (LL) - A (XX) ) *SINCS + A ( L X ) >^ (CCSX 2-S I NX 2) 
A (LL) =Y 
A (XX) =X 

TESTS rCH CCXPL2TI0N 

TEST fCB X = LAST COLUMN 

I? (X-N) 135,140,135 

!1 = X + 1 
GO TO 60 

TEST FOB L = SECCNE xECX LAST CCLUXN 
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r:-(i-(N-i)) Tio, 

L = L+ 1 
GO T C S ? 

If(INC-l) 160,155,160 
IN C=0 
GC TC 50 



CCMFAR2 'IllSESHCLD WITH PINAL NORH 
IF (TliR-AKN.'^X) 165,165,45 

SORT EIGENVALUES ANC EIGENVECTCES 
IC=-N 

BO 185 1=1, N 
IQ=IC+N 

LL=I+ (I*I-I)/2 

JC=N* (1-2) 

DO 185 J=I,N 
JC=JC+N 

NM = J+ (J*J-0) /2 

I? (A (LL)-A (NN) ) 170, 185, 185 

X=A (LL) 

A (LL)=A («K) 

A (MN) =X 

I?(KV-1) 175,185,175 

DO 180 K=1,K 
ILR=IO+K 
IMR=JC+i< 

X= R (ILR) 

R (ILR) =R (lER) 

R (INR) =X 
CONTINUE 
RETURN 
END 



< 
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TYPICAL OUTPUT FORMAT 
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xiESSIAS SAI5IX IS 
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J.[ivj5TLl2b 
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1.0’;51'luO0 




1 . 0C51 tJv3v)0 




2. .)Ci£<;00.') 


HwiST.v-. vat 


3 X 


iU’S, 


-1112.5.432. -- . . 


•4l5.-:>«32 


2 08. -9 697 . 


:.';ai y 


VLCTCr Cf OFTX'IAL X 


OoOOO 


1.JUUC7 1.00017 



siGNr:r- vakiat:on 221 'i = c.iccco 



LCWABL2 IS X VASIID TSDIPISCf STL Y SOS ASSIGSr.L V a I AT J.C ■■. 
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